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INTRODUCTION 

Composites  are  "designer"  materials  in  the  sense  that  a  designer  has  the  freedom  to  prescribe  the 
material  microstructure  such  that  global  response  measures  of  interest  are  optimized  for  a  given 
load  environment  together  with  certain  cost/weight  constraints.  Such  designer  freedom,  however, 
creates  a  need  for  material  response  models  which  are  synthesized  directly  from  material  and 
geometrical  information  at  the  microscructural  level.  This  need,  in  turn,  stems  from  a  desire  to 
avoid  a  myriad  of  experiments  which  may  be  necessary  to  evaluate  the  material  parameters 
associated  with  piicnomenologicai  models. 

A  situation  of  special  interest  to  material  designers  concerns  tiber-reinforeed  polymer  ar.d 
metal-matrix  composites  subject  to  dynamic  load  environments.  Within  the  context  of  such 
materials  and  environments,  response  measures  of  stress  wave  attenuation  and/or  dispersion  are 
often  sought.  For  such  problem  types,  one  of  the  earliest  successful  attempts  to  synthesize  a 
global  response  theory  from  microstructural  information  is  due  to  Achenbach  and  Herrmann 
(1968)  who  formulated  a  higher-order  continuum  model,  known  as  the  "effective  stiffness  theory," 
to  simulate  elastic  wave  motion.  Subsequent  extensions  and  applications  of  this  work  were 
conducted  by  Bartholomew  and  Torvick  (1972),  Hlavacek  (1975),  Achenbach  (1975,  1976),  and 
Aboudi  (1981).  By  modifying  the  original  methodology  A boudi  (1982,  1985)  extended  the  linear 
model  to  account  for  viscoplastic  material  response.  In  parallel  to  the  effective  stiffness  theories, 
attempts  were  made  to  develop  mixture  (multi-phase)  continuum  theories  with  microstructure.  A 
representative  cross  section  of  this  subject  includes  the  works  by  Martin  et  al.  (1971),  Choi  and 
Bedford  (1973),  Hegemier  et  al.  (1^73),  Hegeirier  and  Gunman  (1974),  Nayfeh  (1977), 
Murakami  et  al.  (1979),  Nayfeh  et  al.  (1984),  and  Murakami  and  Hegemier  (1986). 

To-date,  the  foregoing  mixture  theories  have  not  been  extended  to  model  nonlinear  material 
responses  for  arbitrary  wave  motion.  In  view  of  the  potential  modeling  capability  of  the  mixture 
descriptions,  and  in  response  to  a  perceived  need,  a  pf^cedure  is  illustrated  in  this  paper  for  the 
mathematical  consirucnon  of  a  higher-order  mixture  description  of  nonlinear  wave  propagation  in 
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unidirectional  fiber-reinforced  composites.  The  lesuldng  model  incorporates  wave  dispersion, 
wave  anenuation,  localized  plastic  flow,  and  effective  anisotropy.  For  clarity  of  presentation,  the 
example  construction  treated  herein  is  applied  to  a  hexagonal  array  of  fiber  reinforcement  and  rate- 
independent  elastoplastic  material  nonlinearities.  While  the  model  construction  procedure  is 
applicable  to  arbitrary  fiber  layout,  rate-dependent  material  response,  and  interfacial  slip,  extension 
and  investigation  of  such  cases  are  deferred  to  later  publications. 

The  model  construction  method  treated  herein  is  based  upon  a  multi-scale  homogenization 
technique  developed  by  Hegemier  and  Gunman  (1974)  for  waveguide  propagation  and  by 
Murakami  and  Hegemier  (1986)  for  arbitrary  linear  wave  motion.  The  methodology  yields  the 
equations  of  motion,  the  appropriate  initial  and  boundary  conditions,  and  a  set  of  consistent  rate 
constitutive  relations.  The  derived  continuum  mixture  theory  is  nonphenomenological  in  the  sense 
that  the  model  is  synthesized  from  the  composite  "microstructure"  which  consists  of  the  fiber 
material  and  geometrical  properties,  the  interface  properties,  and  the  matrix  propenies. 

A  considerable  effort  is  made  to  validate  the  resulting  model  subsequent  to  its  derivation. 
For  this  purpose,  a  numerical  experiment  is  employed  as  the  "exact"  basis  for  comparison.  Here 
numerical  predictions  from  the  continuum  description  are  compared  with  detailed  finite  element 
(FE)  results  for  several  key  time  dependent  boundary  value  problems.  This  task  necessitated  the 
development  of  a  special  FE  code  for  the  model.  The  FE  code  DYNA2D  (Hallquist,  1982)  was 
employed  to  generate  the  "exact"  reference  data;  for  this  purpose  a  fine  mesh  was  used  to  explicitly 
model  the  composite  microstruemre. 

In  addition  to  information  concerning  accuracy,  the  validation  calculations  reveal  some 
interesting  features  regarding  response  characteristics  associated  with  wave  dispersion,  wave 
attenuation,  localized  plastic  flow,  and  effective  anisotropy.  These  calculations  also  furnish 
enlightening  information  concerning  computational  efficiency. 
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FORMULATION 


Consider  a  domain  V  with  a  uniaxial  periodic  array  of  fibers  embedded  in  the  matrix,  as 
illustrated  in  Fig.  1.  The  position  vector  x  is  described  with  respect  to  a  rectangular  Canesian 
coordinate  system  Xj,  i=l-3  with  Xj  in  the  axial  direction  of  the  fibers.  In  the  X2,  X3-plane,  a 
typical  cell  that  represents  the  geometrical  microstructure  of  the  composite  is  shown  in  Fig.  2  for  a 
hexagonal  array.  For  notacional  convenience  forms  ( a=l,2  denote  quantities  associated  with 
material  a;  here  a=l  represents  fiber  and  a=2  represents  matrix.  The  notation  V  is  the  gradient 
operator  with  respect  to  x  and  (* )  =  a(  )/3T  will  be  employed  in  which  t  represents  time. 
Furthermore,  overbars  will  designate  dimensional  quantifies;  the  lack  of  overbars  will  indicate  non- 
dimensional  quantities. 

The  governing  relations  for  the  displacement  vector  IT  t=u;)  and  the  stress  tensor  a  (=d;j)  in 
each  material  domain  are: 

(a)  Equations  of  motion 

+  =  inV®  (1) 

where  (=fi(“^)  is  a  constant  body  force,  denotes  mass  density,  and  (  )’f  denotes 
transposion  of  the  tensor  ( ); 

(b)  Rate  consfitufive  relations 


•  Tl  1“' 

in  V 


12) 


e  U  =  1  {  V  .  u  ^“^  +  (  V  .  u  }  in  V 


(3) 


where  (=eijl“))  is  the  rate  of  deformation  under  the  small  strain  assumption,  and 
(=Ci^il®P))  is  a  tangent  modulus  tensor  which  becomes  a  constant  tensor  C  for  elastic  response; 


(c)  Interface  continuity  relations 

-d) 


u 


-  u 


=  0  ,  .  ■  y  3  -  u*  )  =■  'j  Oil  A 


(4) 
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where  Ai  is  the  interface  between  fiber  and  matrix,  v(')  (=V|(i))  denotes  the  outward  normal  on  Ai, 
and  0  is  the  zero  vector, 

(d)  Displacement  boundary  data  on  and  traction  data  on  dVj  where  3V=9V^uaVj  is  the 
boundary  of  V; 

(e)  Initial  conditions  at  T  =0. 

In  eqns  (l)-(4),  the  Cartesian  components  of  V.ff,  C;e,  and  V.a  are  and 

respectively.  The  initial  boundary  value  problem  defined  by  the  relations  (a)-(e)  on  V  =  o 
V^“))  is  weU  posed. 

Most  domains  of  practical  interest  contain  a  mulrirude  of  fibers;  for  such  domains,  a  direct 
numerical  finite  element  solution  becomes  intractable  even  with  the  use  of  supercomputers.  In  an 
effort  to  alleviate  this  problem,  a  higher-order  mixture  model  is  developed  to  describe  the  average 
deformation  for  both  fiber  and  matrix  simultaneously,  along  with  higher  order  micro-structural 
deformations.  This  procedure  was  successfully  applied  to  elastic  response  of  fiber-reinforced 
composites  with  a  hexagonal  array  of  fibers  by  Murakami  and  Hegetnier  (1986).  In  what  follows, 
the  above  model  is  extended  to  include  inelastic  response  of  the  constituents,  eqn  (2). 


MODEL  DEVELOPMENT 

The  multivariable  field  representation 

The  derivation  of  the  model  commences  with  a  scaling  of  both  dependent  and  independent 
variables.  To  this  end,  it  will  be  convenient  to  nondimensionalize  the  basic  equations  by  using  the 
following  quantities  (Hegemier  and  Gurtman,  1974): 
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A  : 

typical  macrosignal  wavelength 

A  : 

typical  fiber  spacing  or  cell  dimension 

^(m)  ’  P(m) 

reference  wave  velocity  and  macrodensity 

E  =  p  ■ 

^(m)  F(ni)  '“(m) 

reference  modulus 

t.  ,  =  A/C,  .  : 

(m)  (m) 

typical  macrosignal  travel  time 

£  =  A/ A  : 

ratio  of  micro-to-macrodimensions 

With  the  aid  of  the  above  notation,  nondimensional  variables  are  introduced  according  to: 

(x,  u^“^)  =  (x,  u^“^)/A,  t  =  t/i,^,, 


(C,  (y)  =(C,  (J)  p  =p  /p 


The  periodicity  of  both  fiber  array  structure  and  material  propenies  define  a  cell  in  the 
X2,X3-plane  as  shown  in  Fig.  2.  The  field  variables  in  the  composite  exhibit  significant  variation 
over  two  length  scales:  the  global  and  cell  geometry.  Funher,  an  order  of  magnitude  difference  in 
the  two  length  scales,  suggests  the  use  of  a  multiscale  or  multivariable  asymptotic  technique 
(Babuska,  1976;  Tartar,  1977;  Benssousan  et  al.,  1978;  Hegemier  et  al.,  1979;  Sanchez-Palencia, 
1980;  Murakami  et  al.,  1981).  One  introduces  microposition  vector  x*: 

x*  =  x/e  (6) 


Field  variables  are  now  considered  to  be  functions  of  both  macro-and  microposition  vectors: 

q  (  x,  t )  =  q  (  X,  X*,  t ;  £  )  (7) 

where  xe  V  and  x*e  A(“).  The  cell  domain  is  heterogeneous  in  the  x*-space  and  consists  of  A(') 
and  A(2)  occupied  by  the  fiber  and  matrix,  respectively;  the  macrodomain  V  in  the  x-space  becomes 
homogeneous  and  is  shared  by  the  two  constituents.  Homogeneity  of  geometrical  and  material 
properties  in  the  xi-direction  eliminates  xj*  dependence  in  eqn  (7);  heterogeneity  is  manifested 
only  in  the  X2,  X3-plane.  Consequendy,  spatial  derivatives  take  the  new  form: 
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(8) 


„  (a)  „  -  (a)  ^  1  r7_  ~  (Cl) 

Vq  =Vq  +  —  V^q 

£ 

where  V*  is  the  gradient  operator  with  respect  to  x*,  and  ( ),i»  =3(  )/dxi*=0.  In  the  sequel  q  will 
be  written  as  q  for  notational  simplicity. 

The  operations  (7)  and  (8)  when  applied  to  all  field  variables,  nondimensionalized  by  eqn 
(5)  lead  to  the  following  synthesized  field  equations: 

(a)  Equations  of  motion 

V  .  a  <“>  +  i  V .  ®  f  =  p  u (j'“U  in  V  and  a'“’  (9) 

£ 

(b)  Rate  constitutive  relations 

<j  (“^  =  C  (  e  (  u  +  -!-  e*  (  u  }  in  V  and  (10) 

£ 

where 

e^u^“^)=l{  Vu^“^  +  (Vu^“y  )  (11a) 

e*(u^''^)=i{  V*  u  (  V*  u  )  (11b) 


(c)  Interface  continuity  relations 

(2)  ..  (1) 


u'  =0, 


‘<'’.(s®-a<‘>)  =  0  on  A, 


where  V*(^)  is  a  unit  outward  normal  to 

(d)  Displacement  boundary  data  on  and  traction  data  on  dVj; 


(12) 


(e)  Initial  conditions  at  t=0. 

In  eqns  (9)-<ll)  it  is  understood  that  (  ),i*=0.  The  synthesized  field  variables  (7)  are  now 
continuous  with  respect  to  x  in  V  and  may  be  piecewise  continuous  with  respect  to  x*  in  the  cell 
due  to  the  heterogeneity  of  the  composite. 

At  this  point,  the  variation  of  field  variables  which  satisfy  the  periodicity  with  respect  to  x* 
is  assumed.  According  to  this  condition  field  variables  take  equal  values  on  opposite  sides  of  the 
cell  boundary.  Let  the  fundamental  translation  vectors  of  the  periodic  array  in  the  xi,  X3-plane  be 
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denoted  by  ed^  and  ed2.  In  the  X2*,  X3*-plane  dj  and  d2  become  the  fundamental  translation 
vectors  (for  example,  Kittel,  1971).  Employing  the  direct  notation  x=(xi,  X2,  X3),  x*=(0,  x-.*, 
X3*)  the  x*-perio<licity  condition  for  a  general  cell  with  the  above  fundamental  translation  vectors 
is  expressed  as 

q  (  X,  X*,  t )  =  q  (  X,  X*  +  mj  dj  +  m^  d^,  t )  (13) 


where  m^  and  m2  take  ±1  or  0. 

The  cell  domain  in  the  X2*,  X3*-plane  consists  of  subdomains  A(^)  and  A(-).  Let  the  volume  (area) 
fraction  of  material  a  be  denoted  by  n(“);  it  satisfies 


(1) 


n^^-^  =  1 


(14) 


For  a  hexagonal  array  the  actual  cell  may  be  modeled  as  two  concentric  cylinders  without  loss  of 
accuracy  in  dispersion  spectra.  For  the  concentric  cylinders  model  the  cell  subdomains,  A(L  and 
A(2)  ,  are  represented  as 


=  (  (r,0)  1 

0  <  r  <  >/n^  ,  0  <  0  <  2;t  } 

(15a) 

=  {  (  r ,  0  )  1 

>/n^  <  r  <  1  .  0  <  0  <  2:t  } 

(15b) 

where  (r,  0)  are  polar  coordinates  defined  in  the  X2*,  X3*-plane  such  that 

+  x*^  ,  tan  0  =  x*j  /  x*^ 

(15c) 

For  the  concentric-cylinders  model  the  cell  boundary  is  denoted  by  r=l  and  the  periodicity 
condition  (13)  simplifies  as  follows; 

q(x,r,0,t)  =  q(x,r,0-ni,t)  onr=l  (16) 


When  Fourier  transforms  are  applied  to  both  spatial  variable  x  and  time  t,  the  x*- 
periodicity  (13)  takes  the  same  form  as  the  Floquet  and  Bloch  theorems  for  hanmonic  wave  in 
periodic  structures  (Brillouin,  1946;  Kohn  et  al.,  1972).  Although  equation  (13)  compromises  the 
ability  to  capture  micro-boundary  layer  effects  on  a  cell-scale,  it  provides  an  economical  model  for 
predicting  global  boundary  layer  effects  on  a  scale  of  down  to  a  few  cell  lengths  (Murakami, 
1990),  which  is  sufficient  in  most  problems  of  interest 
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Weighted  residual  procedure 

In  this  subsection,  a  weighted  residual  procedure  is  introduced.  Thi.s  procedure  will  be 
subsequently  used  to  eliminate  x*  from  all  field  variables  through  an  averaging  operation,  and  to 
establish  appropriate  equations  of  motion  for  the  resulting  avenge  fields. 

To  begin,  let  a=l,2,  denote  the  space  of  all  H^-funciions  (for  example,  Hughes, 
1987)  q(x,  X*;  t)  on  V  with  respect  to  x  and  on  A(“)  with  respect  to  x*  that  are  x‘-periodic 
according  to  eqn  (13).  Functions  qlfi  and  ql^l  may  suffer  a  discontinuity  on  the  interface  Ar.  .Any 
vector  ul“)  whose  components  with  ul“)=u(“)  on  where  ul“^  is  the  specified 

boundary  displacement  vector,  will  be  called  an  admissible  trial  displacement.  Any  function 
whose  component  5ujl“)e  with  5u(“)=0  on  will  be  called  a  weighting  function  or  an 


admissible  variation  of  Ui(“). 

Next,  consider  the  weighted  residual  R  defined  by 


jlX  j(V.a<“>  +  iv. 


■ijv 


V  a=l  (a) 
A 


y.®.*<->.5u  <->ds*A 


i  r[iv*<'>  (a'->-a"V(5u"C6uA 

C  J 


+  (  (  5u  -  5u  )]  ds*]  dV 

+  J  (  I  ( ''T  -  V  •  <j  ^“^)  •  5u  dA*  )  dA  =  R  (17) 

3V.^  a=l  (a) 

A 


where  dV=dxidx2dx3,  dA*=dx2*dx3*,  9A  is  the  cell  boundary,  ds*  is  an  infinitesimal  line 
element,  V*(2)  is  a  unit  outward  normal  to  Al^),  and  'T(“)  denotes  the  traction  vector  acting  on  an 
infinitesimal  surface  element  dA  with  a  unit  outward  normal  V  on  DVj.  By  vinue  of  the  x'*- 
periodicity  the  integrations  with  respect  to  macrocoordinates  x  are  carried  out  over  the  entire 
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domain  V,  whUe  that  with  respect  to  the  microcoordinaies  x*  is  performed  over  the  cell  subdomam 
A(“). 

If  R=0  is  sadsfied  for  all  admissible  5u(“)  which  satisfy  (13)  and  are  arbitrary  over  V,  on 
dV-r,  and  on  then  it  is  evident  that  weak  soludon  of  the  local  equations  of  motion  (9)  have 
been  generated.  These  weak  soludons  then  sadsfy  the  traction  boundary  condition  specified  in  (d) 
on  dVj  and  the  traction  continuity  (12)  on  A[. 

From  (17)  with  R=0,  Gauss'  theorem,  and  the  x*-periodicity  condition  (13),  one  obtains 

V  a=l  (a)  ^ 

A 

2 

+  1  J(  5u  -  Su  ^^^)  ■  T*  ds*  ]  dV  =  J*  Jsu  "'T  cLA*  )  iL-\  (IS) 

^  A,  8Vt 

V.  here  the  component  of  5e:ij  is  Se^aij,  and 

5e^“^=l(  V5u^“V(V5u^“y},  =  i  (  V*  5u  (V*  5u  ^“')'^}  (19) 

Equation  (18)  can  be  envisioned  as  the  principle  of  vinual  work  for  the  synthesizeu  fields. 
This  principle  furnishes  a  useful  tool  for  generating  the  equations  of  motion  associated  with  any 
order  of  continuum  models. 

Asymptotic  analysis 

In  Older  to  generate  a  continuum  model  from  (18),  the  assumed  x'^-dependency  of  the 
displacement  field  must  be  described  explicitly.  The  necessary  microscructural  information  for  this 
operation  was  obtained  for  elastic  response  of  composites  with  a  hexagonal  array  of  fibers  by  an 
asymptotic  procedure  (Murakami  and  Hegemier,l986).  This  procedure  is  based  upon  the  premise 
that  the  typical  cell  length  is  much  smaller  than  the  macrodimension,  e«l.  Therefore,  the  form  of 
scaled  eqns  (9>-(12)  suggest  the  expansion  of  the  dependent  variables  in  the  asymptotic  series 
(Lene  and  Leguillon,  1982): 


u  X  ,  XM  .  e  )  =  u  X  ,  X*.  t )  +  E  u  (  X  .  X*.  t ,  ^  e  '  u  (  x  ,  x*.  t  j  -  ...  (20a) 
a  '  ^  (  X ,  X*,  c  ,  c  )  =  -  <y  (  X  ,  x*.  t )  +  O'  (  x  ,  x*.  t  )  -r  e  a  (  x  ,  x*.  t  j  -r  ...  ^20b; 

~^(0)  ■*■^'-(1)  ^  '-CJ  ^  •••  (-0^'-* 

where  U(n^(-),  <J(n)^~^  and  satisfy  the  x*-pericxliciLy  cond'don.  In  the  sequel,  a  class  of 

hardening  elastoplastic  materials,  which  admit  a  ."ate  potential  and  have  positive  definite  tanaent 
modulus  tensor,  is  considered.  Consequently,  C(^0^(‘a)(cp)  assumed  to  be  svmmetr.c  and 

posiuve  definite  in  tlic  expansion  (20c);  this  obviously  holds  for  elastic  response;  where  only 
=  C  (“)  is  required. 

If  eqn  (20)  is  substituted  into  eqns  (9)-(12)  and  the  coefficients  of  different  powers  of  e  are 
equated  to  zero,  a  sequence  of  micro-boundary  value  problerr:  (.\lBVPs)  defined  on  the  cell  is 
obtai — 1.  The  first  three  sets  of  MBVPs  for  the  coefficients  of  E'2,  £-1,  and  E^  are  defined  in  what 
follows. 

ME  VP  for  0(e-2): 


=C 


k  ^  *  A  (Ot) 

'  ^  (.1)  "0  ^  ^ 

t21a) 

.nA- 

(21bj 

A) -""’P""  '■"A, 

(21c, 

From  eqns  (21a,b)  the  operator  for  u,q)(“)  may  be  expressed  as 

^  ^  “  (0)  ^  ^  ^  (0)  •  ^  ^  (0)  ^  J  ® 

A  solution  of  the  problem  is  n^Q)  which  is  independent  of 

^  (0)  “  (0)  ^  ^  e*  (  u  =  o  ^  =  0 


MB  VP  for  0(E'^): 


v*.c;“;=o  in  A 


=  0  on  A, 


(24b) 


(24c.d) 


Equations  (24a,b)  imply  that 


L  (  u =  v.  (  C  e«(  u  “> )  1  =  -  v.  I  C e  (  u  )  (25) 

Equation  (25)  shows  that  is  governed  by  the  same  operator  as  that  of  eqn  (22)  for 

except  for  the  right-hand  side  (RHS).  Even  if  uhe  RHS  of  eqn  (25)  is  nonzero  it  vanishes  when 
integrated  over  the  cell.  As  a  result,  the  integrability  condition  for  is  satisfied.  The  form  of 

the  forcing  term  in  eqn  (25)  suggests  the  following  expression  for  U(,)l“h 

u;“;(x,x^t)  =  ep^(u^O))x”'"^’^^)  ^-6) 


where  xP4  is  x*-periodic.  The  substitution  of  (26)  into  (24)  yields  a  .\IBVP  for  each  xP4;  this  is 
continuous  over  the  cell  due  to  the  perfect  bond  condition  (12).  These  problems  are  defined  up  to  a 
constant  vector  with  respect  to  x*.  This  constant  term  may  be  included  in  0  ■ 

Therefore,  it  is  convenient  to  choose  xP4  such  that  its  integration  over  the  cell  vanishes: 


If 

0=1  (u) 

A 


dA*  =0 


MBVPfor  0(£^): 


-.C;“;^^P^;{e(U(0))  +  e*(  □;“;)}  in  A^“^ 


(2Sb) 


(28c,d) 


At  this  poinq  it  is  instructive  to  outline  0(1)  homogenization  procedures  to  compare  it  with 
the  proposed  0(e)  homogenization  procedure.  Both  homogenization  procedures  require  the 


solution  of  the  MB  VPs  for  defined  by  eqns  (24)  and  (26).  The  0(1)  equations  of  mouon 

are  obtained  by  imposing  the  integrability  condition  for  on  eqn  (28a)  without  solving  the 

MBVPs  for  (Benssousan  et  al.,  1978;  Sanchez-Palencia,  1980).  .According  to  the 

Fredholm  alternative  theorem,  the  problem  defined  in  eqns  (28)  has  a  unique  solution  up  to  a 
constant  vector  with  respect  to  x*,  if  the  operator  for  Uj^2)*-“^  satisfies  the  integrability 

condition —  the  range  of  L(U(2'/“^)  is  orthogonal  to  its  kernel  t)  (for  example, 

Marsden  and  Hughes,  1983).  The  same  0(1)  equations  of  motion  can  be  obtained  by  substituting 
the  trial  displacenient  field  (31)  in  eqn  (18)  and  keeping  only  0(1)  terms.  As  a  result ,  the  0(1) 
model  neglects  the  kinetic  energy  associated  with  the  0(e)  displacement  and  fails  to  model 
harmonic  wave  dispersion.  This  deficiency  is  improved  in  the  0(£)  homogenization  procedure. 


Trial  displacement  and  mixture  equations  of  motion 

The  development  of  an  0(e)  homogenization  commences  with  the  definition  of  an  average 
displacement  field  for  each  constituent,  keeping  terms  up  to  and  including  O(e-)  terms: 

(  X,  t )  =  -L  I  {  U(0)  +  £  u|“J  +  e2  u|“J  )  dA-  (29) 

A  (a) 

A 

where  A(“)  denotes  the  cell  subdomzun  for  integration  and  the  area  for  algebraic  operations. 
Tquatiou  (25)  shows  that  is  excited  by  tip(O)«q'^^q(0)'p-  Therefore,  the  mixture  formulation 

becomes  more  tractable  by  introducing  generalized  displacement  variables  (parameters)  which 
represent  Up(0).q+Uq(Q’.p,  such  that 

S  (X,  t)  =  S  =— f  u^“^V*^^^ds*=i.  fu5“J-V*^^^ds»  (30) 

^  qp  aJ  0) 


where  A=A(ri+A(2)  denotes  the  total  area  of  the  cell. 

This  yields  the  following  trial  displacement  field 

u  X,  xM;  e  )  =  U  x,t)  +  eSp^(x,t)X^‘’(x*)  (31) 

In  eqn  (31)  is  the  average  displacement  associated  with  each  constituent,  while  Spq(x,  t)  is  the 


generalized  displacement  which  represents  the  amplitude  of  the  0(£)  displacement  microstmcture. 
In  the  sequel,  equation  (31)  will  be  used  to  obtain  mixture  equations  of  motion  from  eqn  (18). 


In  order  to  find  the  0(£)  displacement  microstructure  xW(x*)  one  must  solve  six  MB  VPs 
defined  by  eqns  (24)  and  (26).  These  problems  were  solved  analytically  by  Murakami  and 
Hegemier  (1986)  for  a  hexagonal  cell,  consisting  of  elastic  constituents,  approximated  by  the 
concentric-cylinders  model.  The  exact  solution  indicates  a  good  approximation  for  the  0(£) 
displacements,  and  the  following  trial  displacement  field  (in  component  form)  was  constructed  for 
hexagonal  cells: 

u  X  ,  X*,  t ,  £  )  =  U  (  X  ,  t )-!-  £[  S  .  2  (  X  .  t )  cos  0 +  S  j  ^  .  t )  sin  9  ;  g"^\r)  (32a) 


where 


g<»'(r)  =  (-ir'  JL(r-5  i) 


(a) 


az  r 


S  =  S 
■^23  '^3  2 


(32b) 

(32c) 


and  where  is  the  Kronecker  delta.  The  generalized  displacements  Sjj  are  not  displacement 
compoments  but  parameters;  it  is  convenient  to  employ  the  component  form  in  eqns  (32).  The 
functions  g(“)(r)  cos9  and  g(^)(r)  sin0  are  the  approximations  for  and  satisfy  the  x*- 
periodicity  condition  (16)  and  the  normalization  condition  (27).  The  effectiveness  of  the  above  trial 
displacements  to  simulate  harmonic  v/ave  dispersion  was  demonstrated  for  hexagonal  and  square 
arrays  in  the  above  reference.  For  arbitrary  cells  and  elastic  constituents  one  can  num.erically  solve 
the  MBVPs  for  XP9  by  finite  element  methods  and  numerically  construct  approxim.ate  solutions. 
These  approximate  solutions  are  functions  ofx*  and  independent  of  material  properties;  therefore, 
they  apply  to  both  isotropic  and  orthotropic  constituents. 

For  nonlinear  responses  jN  can  be  found  by  solving  the  rate-MBVPs  since  the  tangent 
moduli  must  be  evaluated  for  each  e(U(0));  dus  implies  that  xW  differs  for  each  load  increment.  To 
render  the  following  analysis  tractable,  an  approximate  solution  for  which  is  independent  of  the 
load  increment  was  constructed.  By  virtue  of  very  weak  anisotropy  introduced  to  the  nonzero  and 


off-diagonal  entries  of  elastic  modulus  tensor,  it  is  found  that  the  approximate  foe  elastic 
response  furnishes  a  good  approximation  even  for  elastoplastic  deformation.  This  situation  is 
similar  to  the  nonlinear  plate  and  shell  analyses  in  which  linear  variation  of  in-plane  displacements 
over  the  thickness  of  the  plate  and  shell  -  found  for  elastic  responses  -  yields  a  good 
approximation  even  for  nonlinear  response.  The  soundness  of  the  above  approximation  will  be 
examined  in  the  sequel  by  comparing  the  model  prediction  with  the  one  of  detailed  finite  element 
analyses.  In  what  follows,  the  above  theory  will  be  applied  to  a  hexagonal  cell  with  the  concentric- 
cylinders  approximation. 

Substitution  of  eqns  (32)  into  eqn  (18)  yields  the  mixture  equations  of  motion  and 
associated  boundary  conditions  together  with  the  inherited  initial  conditions;  they  are  given  in 
component  form: 


(a)  Mixture  equations  of  motion 
_  (ct)  _  (cu)  , 


(a). (a)  (a)  (a)  ••  (o) 


n  +(-!)“"*  P  .  n  V  ’  =  n  p  U 

j  I .  j  I  I  ‘ 

e 


M  . ,  ,  -t-  -L  (  a  -  a  -I-  R  T  )  =  I  S  ,  ,  i  =  1,  3 
e 


(34b) 


e 


where  the  average  operations  are  defined  by 


a|““^x,t)=  —  f  a|“^x.x*.  t)dA^ 
‘J  .  (a)  J  ‘J 


*  (a)  J 

A  (a) 

A 


(37) 


2 

e(\‘l.  .,M.j)=  ^  X  J o|“^g^“^cos0,sin0)dA* 


a=l  (a) 

A 


R  , .  =  — i—  f  —  (  a  ?!  cos  20  +  a  sin  20  )  dA*,  i  =  1.  2 
2  1  (2)  .  J  2  '  2 1  3  1 


(38a) 


n  '  'A  r 


R,.= 


1  r  1  ,  _(2)  ,  ^(2) 


COS 20  +  a'*'  sin20)ciA*.  i=l,3  (38b) 


R 


*  afi)sin26dA. 

n^^’A  2r 


3  3 


(38c) 


=I 


a.i  4n'-’ 


(39) 


In  eqns  (36)-(38)  A(=7t)  denotes  the  area  of  the  cell; 
(b)  Boundary  conditions 


U<“>  or 
\ 

n 

specified  for  i  =  1- 3 

(90) 

S.2 

or 

2 

M  j .  V  j  specified  for  i  =  1 , 2 

(41a) 

Si3 

or 

3 

M  j .  V  j  specified  for  i  =  1 ,  3 

(41b) 

^  2  3 

or 

3  2 

(Mj2’+'M.2)Vj  specified 

(41c) 

(c)  Initial  conditions 

1  ’  1 

S  . 

1 

2’^i3’  ^i2’  ^i3  specified  at  t  =  0 

(42) 

It  is  noted  here  that  the  above  mixture  equations  are  identical  to  those  for  elastic  constituents 
(Murakami  and  Hegemier,  1986). 


Incremental  constitutive  relations  and  trial  transverse  stress-rate 

At  fu’st  glance,  it  would  appear  that  the  trial  displacements  (32)  could  be  used  together  with 
the  original  three-dimensional  material  constitutive  relations  (10),  and  the  stress-type  averages 

1  6 


(35)-(38)  to  establish  a  set  of  constitutive  equations  for  the  stress  averages  to  accompany  the 
equations  of  motion  (33)  and  (34).  A  closer  examination,  however,  reveals  that  such  an  approach 
will  not  yield  a  relation  for  the  interaction  body  force  P  in  (36)  and  will  lead  to  a  model  w  hich  is 
too  "stiff'  and  which  exhibits  erroneous  dispersive  characteristics.  These  problems  -  common  to 
the  use  of  direct  variational  methods  -  can  be  alleviated  by  the  use  of  a  judicious  mixed  w  eighted 
residual  procedure  wherein  ihe  trial  functions  include  certain  stress-rate  components  as  well  as  the 
velocity  components.  This  is  an  incremental  version  of  Reissnefs  mixed  variational  principle 
(Reissner,  1984,  1986)  which  was  employed  to  derive  the  mixture  model  for  elastic  constituents 
(Murakami  and  Hegemier,  1986). 

In  order  to  use  the  mixed  weighted  residual  procedure,  it  is  necessary  to  rewrite  the  rate 
constitutive  relations  (10)  in  terms  of  in-plane  strain  rates  and  transverse  stress  rates;  these  are 
shown  in  matrix  form  for  easy  finite  element  implementation: 


«ll=E,.e,,  +  (E,,]{a.}  (43a) 

(ij  =-[E,/c,  i  +  (  o.)  (-iSb) 


where 


fa|=[a  a  a  a  o  1^, 

*■  i-*  22  33  23  31  12-' 

f  e  1  =1  e  e  Oe  2.e  -i-  —  f  e*  e*  1 

jJ  L  ^22  j3  ^^23  ^^31  ^^12^  ~  ^  ^  22  ^33  23  3l  12-' 


In  eqns  (43)  and  (44)  [  is  the  transposition  of  [  ];  subscript  t  denotes  transverse  quantities.  The 
transverse  stresses  are  those  which  appear  in  the  traction  continuity  condition  (12b),  i.e.,  all  stress 
components  except  The  matrices  [Ey]  are  functions  of  the  elements  of  and 

[E22^“^J  is  symmetric  and  positive  definite  for  hardening  materials.  Specific  forms  of  [Eyt^)] 
employed  for  the  numerical  study  are  given  in  Appendix  A  and  obtained  for  the  von  Mises  yield 
criterion  and  associative  flow  rule  with  isotropic  strain  hardening. 

Let  (a=l,2)  denote  the  space  of  all  H^functions  q(x,  x*,  t)  on  V  with  respect  to  x 
and  on  A(“)  with  respect  to  x*  that  are  x*-periodic  according  to  eqn  (13).  Functions  q(*^  and  q^^^ 


may  suffer  a  discontinuity  on  the  interface  Aj.  Any  vector  whose  components  belong  to 
with  on  3V^(“),  where  u(“)  is  the  specified  boundary  velocity,  will  be  called  an 

admissible  trial  velocity.  Any  function  whose  components  belong  to  with  5ui'^)=0  on 
will  be  called  a  weighting  funcdon.  The  space  of  admissible  transverse  stress-rate 
(a=l,2)  consists  of  all  (=1-2)  functions  q(x,  x*.  t)  on  V  with  respect  to  x  and  on  A(“)  with 
respect  to  x*  that  are  x*-periodic.  The  mixed  weighted  residual  procedure  applied  to  the  rate 
boundary  value  problem  defined  by  (9>-(12)  yields  in  matrix  form: 

I[l;i(S'ri°rM5'rr{°rK{5u<“>r{pa<“'i 

V  a=l  (a) 

A 

+  (e;“^}  }  dA^ 

+  1  f[(  (5u^^^}^-(5u^^^}^  {t*}+  {5T*)\{u^‘^)-[u‘^^})]ds-]dV 

=  jl£j(5u  dA*ldA  («) 

9V.J-  CX=1  (aj 
*  A 

where(5u),  (pa),  {T*},  and  ('T)  are,  respectively,  the  matrix  representation  of  5u,  pli,  T*,  and 


For  arbitrary  variation  of  (u(“)}  and  (dj(®)},  one  obtains  the  rate  constitutive  relation  for  {e^(“)}  in 
(43b),  as  well  as  the  rate  equations  of  motion  (9),  the  rate  boundary  conditions  (d),  and  the  rate 
form  of  eqn  (12).  Equation  (43a)  is  considered  to  be  the  definitions  of  The  mixed 

weighted  residual  equation  (45)  with  appropriate  trial  functions  for  {u(“l)  (=u<^“l)  and 
yields  the  rate  constitutive  relations  for  the  stress  averages  in  eqns  (35)-(38).  The  rate  form  of 
(32a)  furnish  a  trail  velocity  field.  The  trail  transverse  stress  rate  has  the  form 

{aj“J,j}(x.x*,t)-t-e  {  a  (  x,  x*,  t )  (46) 


In  accordance  with  the  0(1)  homogenization  procedure,  the  0(1)  transverse  stress-rate  {dt(0)) 
be  constructed  by  using  the  approximate  velocity  field  defined  from  eqn  (31).  For  a  hexagonal 


array,  substituting  the  trial  velcx;ity  field  obtained  from  (32)  into  (10)  an  appropriate  form  is 
obtained  for  the  transverse  stress-rate: 

+  S  1(T(0)]  {t'j(x,t)  (47a) 


a2  2 
r 


where  and  {t  J  are  stress  variables  defined  as 

{T<“>i(x,.)=[,f'  xf;  x<“>f 

T 

{tj}(x,  t)=ft22  ^33  ^23  ^31  ^12^ 


(47b) 

(47c) 


and  where  [T(0)]  is  a  5  x  5  matrix  whose  nonzero  elements  are 

'rM=T,,  =  T2,=-T,,^  =  T,^=T,3  =  cos20 

Tj3=T23=T  3  2  =  T,5  =  -T  5,  =  sin  20  (47d) 

The  0(e)  term  is  governed  by  the  MB  VP  (28)  which  requires  the  solution  of 

The  exact  analysis  of  based  upon  the  expansion  (20a)  and  eqns  (23)  and  (26) 

yields  42  sub-MBVPs.  The  formidable  task  of  solving  for  the  elements  of  U(2)^“^  can  be  alleviated 
by  constructing  approximately  from  (28a).  Examination  of  the  analytical  solutions  for  a 

concentric-cylinders  cell  revealed  that  {cri(i)^“M  which  satisfies 

V*.aJ;;  =  (-l)“"'P/n<“^  inA^“^  (48) 

yields  a  simple  approximation  for  a  general  cell. 

Equation  (48)  is  obtained  by  applying  the  Gauss  theorem  to  eqn  (36)  and  satisfies  the  integrability 
condition  for  U(2)^“^  through  the  explicit  introduction  of  P.  For  a  hexagonal  cell  the  following 
approximation  (adopted  for  elastic  constituents)  is  employed. 

-■ig^“\r)[Q(0)]  {p)(x.t)  (49a) 


where  (P)  is  a  matrix  representation  of  P,  and  [Q(9)]  is  a  5  x  3  matrix  whose  nonzero  elements 


are 


(49b) 


Q22~Q33"3Qi2  2^^*  COS0 

^13“Q32"3'^23“‘2^41~^^'^® 

The  trial  transverse  stress-rate  is  now  obtained  by  substituting  eqns  (47)  and  (49)  into  eqn  (46). 

On  substituting  the  rate  form  of  eqn  (32)  and  eqn  (46)  into  eqn  (45),  one  obtains  the  rate 
form  of  (33),  (34),  (40)  and  (41).  In  addition,  the  arbitrary  variation  of  yields  the  rate 

constitutive  relations  for  the  transverse  stress  variables: 

+  4[ti(i  j  +  ig®’(Ql(p))dA* 


=  r.‘“>  It  (  ( e  J  +  (- 1 _L  {  S ) )  ^  J  (  E 


1  ■=  1  1 


f  4[TftE®l(  g‘->(Q]  {p))d.A- 

a  r  r  ^ 

A 

=  7c[W]{s}  +  J  l[Tf  [Ef^^ef-JdA* 

A 

f  X  J  g^“^[Qr[E^“^]((T^“^}+5^  l[T]{t  j  +  Jg^“^[Q][p))dA^ 


0=1  (a) 

A 


=  7t(  {u^^^})-h7t^([R‘]{s}  j+[r"i(s)  ,+  [  R°Ms}  3) 


0=1  (a) 

A 


where  {  U(“))  is  a  matrix  representation  of  U<“),  h=h(^^-»-  h(-),  and 

[  e  ]  =  r  U  U  U  -i-  U  U  +  U  U  +  U  1 

l*^(U)tJ  ‘■'-^2.2  ^3.3  ^  2  ,3  ^  3 .2  ^3, 1^1, 3  '^1.2  2.  M 

{s}  =  [S22  S33  2S23  Sj3  Sjj] 
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In  eqns  (50)-(52),  [W]  is  a  5  x  5  matrix,  and  [R*],  [R^^],  [R^^^]  are  3x5  matrices.  Nonzero 
elements  of  those  matrices  are 


W,.=-W,,  =  W3  3  =  .^,  W3.=  W33=-W,,=  W3  3=-L 

J.  n  n 


R‘..=RU=R?5=Rm=2,R",=r“=3,  r“3  =  r”3=R3“=r”  =1  ,34) 


The  solution  of  (51)  and  (52)  yields  {t  {t  and  (PJ  in  terms  of  {e(U)t(“)),  (S),  and 

1).  For  inelastic  response  the  integrals  in  (50)— (52)  have  to  be  evaluated  numericaJIv  at 
each  increment  (time  step);  for  elastic  constituents  the  above  relations  can  be  evaluated  explicidy 
and  this  produces  the  results  obtained  from  a  fully  elastic  approach  to  the  problem  (Murakami  and 
Hegemier,  1986).  Substituting  eqn  (47)  into  eqns  (37)  and  (38)  one  fmds 

2  23  3 

M,2  =  3hP2/4,  M33=M23=hP,/4  ,  M33=3hP3/4, 

Mi2=M3j=hP^/2  ,  M22  =  ^23  =  hP3/4,  M3  j=Mj3  =  0  (55) 

Rj,=t,,/n“*,  Rj3  =  (t  33/2  +  t  33)/n">  ,  R3,  =  .i  3,,n'"  , 


The  remaining  constitutive  equations  for  and  Mj  [  are  obtained  from  eqns  (43),  (35)  and 

(37).  The  results  are 

n<“V  a  =  J[  E  e  4  [  E  '“>]  (  (t)  +5^  4  [  t  ]  { i  j  +  i  g'“>[  Q  J  (  P  ))  dA*  (57) 


e’'(Md  =  Z  J«'“’{n}[E‘“;e 


(a) 

1  I 


a=l^(a) 


where 

{^^p}=[Mii  M  j  j  ]^,  {n)=[cos6  sin  9 


(58b) 
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The  above  operations  were  carried  out  at  each  integration  point  in  a  constitutive  subroutine  for  the 
mixture  finite  element  code:  HFEC2D  (Impelluso,  1990). 


MODEL  VALIDATION  STUDIES 

In  this  section,  validation  studies  are  conducted  in  an  effon  to  ascertain  the  simulation 
capability  of  the  0(e)  mixture  continuum  model.  The  model  consists  of  the  equations  of  motion 
(33)  and  (34),  the  boundary  conditions  (40)  and  (41),  the  initial  conditions  (42),  and  the  rate- 
constitutive  relations  (50)-(52)  and  (55)-(58).  The  problems  examined  include  linear  and 
nonlinear  wave  propagation.  In  the  dynamic  response,  mixture  model  predictions  were  compared 
with  experimental  data  or  "exact"  numerical  data  generated  from  DYNA2D  based  upon  detailed 
explicit  modeling  of  fibers  and  matrix.  Material  propenies  of  the  investigated  composite  which 
consists  of  elastic  fiber  and  elastoplastic  matrix  are  shown  in  Table  1. 

For  elastic  harmonic  wave  propagation,  the  model  was  validated  by  comparing  the 
predicted  phase  velocity  spectra  with  experimental  data  (Murakami  and  Hegemier,  1986). 
Therefore,  the  validation  of  the  model  was  conducted  in  the  lime  domain.  As  was  noted 
previously,  the  validation  strategy  is  to  compare  mixture  model  predictions  with  experimental  data 
or  "exact"  numerical  data  generated  from  DYNA2D  based  upon  detailed  explicit  modeling  of  fibers 
and  matrix.  For  this  purpose,  an  explicit  finite  element  code:  HFEC2D  (Impelluso,  1990)  was 
developed  using  four-node  quadrilateral  elements  for  the  generalized  plane  strain  in  the  xi,  xo- 
plane.  The  mixture  element  has  six  nodal  degrees-of-freedom  for  Ui(“)  and  Sj-  (i=I,2).  The 
element  carries  the  microstructure  of  the  cell  at  each  integration  point  where  the  numerical 
integration  of  incremental  constitutive  equations  (50)-(58)  are  conducted.  For  simplicity  of 
notation  in  the  numerical  results,  dimensional  quantities  are  referred  without  overbars. 

The  geometry  of  the  wave-reflect  problem  is  shown  in  Figs.  3  with  meshes  for  HFEC2D 
(Fig.  3b)  and  the  detailed  DYNA2D  (Fig,  3c).  A  composite  half  space  with  free  boundary  at  X2=0 
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was  loaded  uniformly  with  respect  to  X3  under  plane  strain  condition  in  the  xi -direction.  For  this 
globally  one  dimensional  wave  phenomena  in  the  X2-direcrion,  a  column  of  cells  of  width 
A=0.0975  cm,  shown  in  Fig.  3a,  is  discretized  by  the  mesh  shown  in  Fig.  3c  for  DYNA2D 
calculation;  for  the  mixture  model  only  one  row  of  elements  shown  in  Fig.  3b  is  employed.  The 
following  bondary  conditions  were  pwsed  for  the  DYNA2D  calculation; 

a22  =  cyo  {  H(t)-H(t-tQ)  ),  0^3  =  0  atx^  =  0 

03  =  023=0  atX3  =  0,  A  (59) 

For  the  mixture  model,  the  corresponding  boundary  data  was  specified  as 

a^“^  =  n^“^Oo{  H(t)-H(t-g).  o  =  M  ,  ,  =  ,  2  =  0  at  x  ,  =  0 

U  f  ^  =  O  =  M  ^  ^  j  2  =  0  at  X  j  =  0,  A  (60) 

where  H(t)  denotes  the  Heaviside  step  function,  and  to=  3}J.sec.  is  the  pulse  duration. 

A  load  of  Co  =  1  X  10^  dyn/cm^  is  applied  to  induce  a  purely  elastic  response  in  both  constituents, 
while  a  load  of  oq  =  3  x  10^  dyn/cm^  is  applied  to  induce  an  elastoplastic  response  in  the  matrix. 
The  numerical  results  are  shown  for  observation  points  located  in  the  33rd  cell.  The  time 
variations  of  fiber  particle  velocity,  at  r  =  0,  and  of  matrix  panicle  velocity  at  r  =  1  and 
0=0*  are  shown  ,  respectively  in  Figs.  4a  and  4b.  The  corresponding  time  variations  for  the 
elastoplastic  case  are  shown  in  Figs.  5a  and  5b.  Arrival  time,  peak  response,  and  damping  are 
well  correlated  by  the  mixture  elements.  The  dispersive  behavior  is  evident  and  well  matched; 
funhermore,  the  spreading  of  the  wave  pulse  to  a  duration  larger  than  3psec.  is  demonstrated  in 
both  DYNA2D  and  P{FEC2D. 

Special  attention  is  paid  now  to  the  localized  plastic  deformation  in  the  composite.  Figure  6 
shows  the  effective  plastic  strain  contour  obtained  from  DYNA2D;  plastic  deformation  is  localized 
near  0=0*  and  180*.  To  examine  the  capability  of  the  model  in  predicting  those  localized  effects, 
the  effective  stresses  at  0=150*  and  180*  both  at  r=0.761  are  shown,  respectively,  in  Figs.  7a  and 
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7b.  From  those  figures,  the  localization  is  modeled  accurately  by  the  mixture  element.  This 
accuracy  is  also  reflected  in  the  predictions  of  the  effective  plastic  strain. 

In  order  to  further  assess  the  accuracy  of  the  model,  the  waveguide  problem  illustrated  m 
Fig.  8  is  conside.^'ed.  The  load  is  applied  at  the  boundary  xi=0  in  the  fiber  axis  direction. 
Hexagonal  symmetry,  approximated  by  axisymmetry,  allow's  for  an  extraction  of  an  a  cell  of  radius 
A=0.0975  cm.  First,  a  comparison  is  made  witii  a  shock  tube  test  conducted  by  the  Ac’-ospace 
Corporation  and  reported  by  Hegemier  et  al.  (1972).  The  composite  is  subjected  to  a  step  pressure 
loading  of  4.826  x  10^  Pa.  The  results  of  the  model  are  obtained  by  using  1 10  elements.  The 
comparison  of  the  rear  surface  velocities  on  a  specimen  6.33  mm  thick  is  shown  in  Fig.  9.  In  the 
above  expenment  the  nonlinear  effect  was  negligible. 

Numerical  experiments  were  conducted  to  test  the  model's  capability  for  predicting 
elastoplastic  wave  propagation.  The  following  boundary  conditions  with  respect  to  the  cylindrical 
coordinate  system  were  pciCd  for  the  DYNA2D  calculation;  \\  and  x:  are,  respectively  the  axial 
and  radial  coordinates: 

a  ^  j  =  Oq  (  H(t)  -  H  ( t  *  t  q)  ),  a  j  2  =  0  at  x  ^  =  0 

u  2  =  a  j  2  =  0  at  X  ,  =  0,  A  (61) 

For  the  mixture  model,  the  corresponding  boundary  data  was  specified  as 

^(cu)^n<a)^^{  H(t)-H(t-g),  =  M  ,  =,<1  ^  2  =  °  at  x  ,  =  0 

U  2“'  =  ai“^  =  M  1  2  =M2  2  =  ^  at  x,=  0,0.51  A  (62) 

A  pulse  load  oo=  5.0  x  10^  dyn/cm^  of  duration  io=  Spsec.  was  applied.  Axial  velocities  at 
xi=43.5A  were  plotted  in  Fig.  10a  for  fiber  and  in  Fig.IOb  for  matrix. 

The  above  comparisons  with  the  detailed  FE  analyses  indicate  the  cost  efficiency  of  the 
mixture  element  due  to  the  coarseness  of  the  mesh;  the  mixture  code,  HFEC2D,  runs  at  least  t  .i 
order  of  magnitude  faster  than  the  detailed  finite  element  computation. 


CONCLUDING  REMARKS 


The  construction  of  a  higher-order  mixtuie  description  of  fiber-reinforced  composites  has 
been  demonstrated  herein  for  the  case  of  material  nonlinearities.  For  simplicity  of  presentation, 
composites  with  a  hexagonal  array  of  fibers  and  elastoplastic  matrix  were  considered.  The 
methodology  is  based  upon  an  asymptotic  homogenization  method  and  yields  me  equations  of 
modon,  the  appropriate  initial  and  bouiidary  conditions,  and  a  set  of  consistent  rate-constitutive 
relations.  For  transient  response  a  finite  element  wave  code  was  developed  for  the  mixture  model 
to  so've  linear  and  nonlinear  problems;  results  using  this  code  were  compared  with  those  from 
DYNA2D  in  which  a  fine  mesh  was  utilized  to  explicitly  model  the  microstracture  ot  the 
composite.  These  comparisons  reveo’  that  the  mixture  model  is  capable  of  furnishing  an  accurate 
and  economical  description  of  complex  wave  phenomena.  In  the  foregoing  analyses  the 
imponance  of  wave  dispersion  and  attenuation  effects  was  confirmed  for  nonlinear  as  well  as  linear 
composite  responses. 
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-  APPENDIX.  DEFINITIONS  OF  [Ejj]  IN  EQN  (43) 

-  -  “It  is  computationally  advantageous  to  revmte  eqn  (10)  in  terms  of  the  elastoplastic  comptiance 
matrix  D=C*^  For  a  von  Mises  yield  criterion  and  associated  flow  rule  with  linear  strain 
—hardening,  D(“Kep)  may  be  expressed  as  - 


(a)(ep)  _  p  (a) 


4  a  H'^“^ 


^(a)  ^(a)T 


(Al) 


In  eqn  (Al),  H'  is  the  strain  hardening  parameter,  s  (=a-tra.5)  is  the  deviatoric  stresses,  and  a 
(=V3s:si/2)  is  the  Mises  effective  stress. 


Rewriting  the  rate  compliance  relation  by  using  a  6  x  6  matrix  [D]  for  D,  one  finds 

(a) 

1  1 


- ^13  ^14  ^15  ^16^ 


D 


(a)(ep) 
1  1 


D 


(a)(ep) 
1  I 


i,j=2.6 

^11 


(A2) 
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Table  1  Material  Properties  for  Wavereflect  Problem. 


Material 

Volume 

Fraction 

Density 

Young's 

Modulus 

Poisson's 

Ratio 

Yield 

Stress 

Hardening 

Parameter 

a 

n(“) 

p(“)(g/cc) 

E(“Hdyn/cm“) 

v(«) 

aY(dyn/cm‘) 

H'(dynycm-) 

1 

0.272 

1.85 

292.0  (109) 

0.3776 

. 

2 

0.728 

1.29 

82.24  (109) 

0.35*? 

1.33  (109) 

11.38  (109) 
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reinforced  composite  domain  V  with  a  uniaxial  periodic  array  of  fibers. 
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Fig.  4.  Time  variation  of  particle  velocity  for  the  elastic  wave-reflect  analysis:  (a)  fiber  -  U2ri) 
(b)  matrix  - 
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Fig.  5.  Time  variation  of  particle  velocity  for  the  elastoplastic  wave-refiect  analysis:  (a)  fiber  -  U2ri) 
(b)  matrix  - 
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Fig.  6.  Effective  plastic  strain  contour  obtained  from  DYNA2D. 
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Fig.  7.  Time  variation  of  effective  stress  in  the  matrix  domain  for  the  elastoplastic  wave-reflect 
analysis:  (a)  9  =  150*,  r  =  0.761,  (b)  0  =  180*,  r  =  0.761. 


